
############################################
# Replication Code                         #
# ISQ Accepted Manuscript: 2016-02-0120-R2 #
# Prepared by Sarah Dreier                 #
# September 2017                           #
############################################

rm(list=ls())
library(Hmisc)
library(dplyr)

########################################
########################################
# Figures 1 A-C Pentecostal Size #######
# Data: World Religion Project,  #######
#       World Christian Database #######
########################################
########################################

#World Religion Project (Correlates of War): National Data
wrp <- read.csv("wrp_national_africa.csv", header=TRUE, sep=",", stringsAsFactors=FALSE)

#Subset to Sub-Saharan Africa
wrp <- subset(wrp, Country!="Morocco" & Country!="Algeria" &
                    Country!="Tunisia" & Country!="Libya" &
                    Country!="Egypt" ) #keep Sudan

#Subset to 1974 - 2010
wrp <- subset(wrp, Year>1974)

#Calculate yearly statistics
year_data <- split(wrp, wrp$Year)

year_stats_fct <- function(dat){
  
  yr_protestant_pct <- mean(dat$chrstprotpct)
  yr_protestant_total <- sum(dat$chrstprot) 
  
  yr_ang_pct <- mean(dat$chrstangpct)
  yr_ang_total <- sum(dat$chrstang)
  
  yr_all_christ_pop <- sum(dat$Christ.Total)
  yr <- dat$Year[1]
  
  return(  data.frame( yr_protestant_pct, yr_protestant_total,
                       yr_ang_pct, yr_ang_total, 
                       yr_all_christ_pop, yr)   )  
}

agg_data <- lapply(year_data, year_stats_fct)
agg_data <- bind_rows(agg_data)

agg_data$christ_perc_prot <- agg_data$yr_protestant_total / agg_data$yr_all_christ_pop
agg_data$christ_perc_ang <- agg_data$yr_ang_total / agg_data$yr_all_christ_pop

# World Christian Database (World Christian Research)
wcd <- read.csv("wcd_ssa.csv", header=TRUE, sep=",", stringsAsFactors=FALSE)

# Rename variables: all numbebers in columns 6-14 are percentage of overall Christian population
names(wcd) <- c("ctry", "yr", "pop", "chr_pop", "perc_chr", "chr_growth", "orth",
                "cath", "ang", "prot", "indep", "unafil", "double_affil", "evan", "renewalist")

# Subset to SS Africa Continent-wide data
wcd_sub <- data.frame(wcd[4:5,])

######################################################
# Figure 1A: Pentecostal Size (% overall Population) #
# Data: World Religion Project                       #  
######################################################

pdf(file="Fig_1A_Protestant_Size1.pdf",width=5,height=4)
with (agg_data, plot(yr, yr_protestant_pct, ylim=c(0,.15), type="l", lwd="2", bty="n",
                     xlab="", 
                     ylab="Percent of SS African population",
                     main="Protestant and Anglican size,\nrelative to total population\n(Sub-Saharan Africa)", 
                     cex.main=.9) )
with (agg_data, lines(yr, yr_ang_pct, lwd="2", col="azure4") )
text(2005, .113, "Protestant\n(including Evangelicals)", cex=.8)
text(2005, .03, "Anglican", cex=.8)
mtext("Source: World Religion Project (Correlates of War)", side=1, line=2.5, cex=.75)
dev.off()

###############################################################
# Figure 1B: Pentecostal Size (compared to other Protestants) #
# Data: World Christian Database                              #
###############################################################

pdf(file="Fig_1B_Protestant_Size3.pdf",width=5,height=4)
with(wcd_sub, plot(yr, ang, xaxt="n", xlim=c(1945, 2040), ylim=c(0,40), bty="n", pch=16, 
     type="b", lwd="2", col="azure4",
     ylab="Percent of SS African Christian population", xlab="",
     main="Protestant, Anglican and Pentecostal size,\nrelative to Christian population\n(Sub-Saharan Africa)",
     cex.main=.9) )
text(2025, 10, "Anglican", cex=.8, col="azure4")

points(wcd_sub$yr, wcd_sub$prot, col="azure4", type="b", pch=16, lwd="2")
text(2026.5, 28.5, "Protestant\n(mainline)", col="azure4", cex=.8)

points(wcd_sub$yr, wcd_sub$renewalist, type="b", pch=16, lwd="2")
text(2027.3, 35.7, "Pentecostal\n(Renewalist)", cex=.8)

axis(1, at = c(1945, 1970, 2015, 2040), labels = c("", "1970", "2015", ""))
mtext("Source: World Christian Database", side=1, line=2.5, cex=.75)

dev.off()

################################################################
# Figure 1C: Pentecostal Size (% overall Christian Population) #
# Data: World Religion Project                                 #
################################################################

## Not used in final manuscript

# pdf(file="Fig_1C_Protestant_Size2.pdf",width=5,height=4)
# with (agg_data, plot(yr, christ_perc_prot, ylim=c(0,.35), type="l", lwd="2", bty="n",
#                      xlab="", 
#                      ylab="Percent of SS African Christian population",
#                      main="Protestant and Anglican size,\nrelative to Christian population\n(Sub-Saharan Africa)",
#                      cex.main=.9) )
# with (agg_data, lines(yr, christ_perc_ang, lwd="2", col="azure4")  )
# text(2005, .255, "Protestant\n(including Evangelicals)", cex=.8)
# text(2005, .105, "Anglican", cex=.8)
# mtext("Source: World Religion Project (Correlates of War)", side=1, line=2.5, cex=.75)
# dev.off()

##########################################
##########################################
# Figures 2 & 7 from Anglican Data       #
# Data: Compiled by S. Dreier, 2013-2017 #
##########################################
##########################################

ang_base <- read.csv("ang_data.csv", header=TRUE, stringsAsFactors=F)

ang <- subset(ang_base, total_church_body=="total") #subset to church body, not country

############################################################
# Figure 2: Anglican Membership Size and Relationship w US #
############################################################

ang <- ang[order(-ang$church_pop_2016),]
ang$ctry_pop <- ang$country_pop_2016 / 1000000
ang$church_pop <- ang$church_pop_2016 / 1000000
x <- c(ang$names[1:7], rep("", 5) )

pdf(file="Fig_2_Ang_Pop_Relationship.pdf",width=5,height=4)
par(mai=c(.9,1.1,1,.5))
plot(ang$church_pop, ang$us_relations, cex.main = .9, ylim = c(.8, 3.3), xlim = c(0, 20), bty="n", 
     xlab = "Church population (millions, 2016)", ylab = "", yaxt="n",
     main = "African Anglican Church relations with US Church\n(by church population)", cex.main=.95
)
axis(2, at = 1:3, labels = c("Relationship\nMaintained", 
                             "Relationship\nStrained",
                             "Relationship\nEnded"), las=2, cex.axis=.7)
abline(lm(ang$us_relations ~ ang$church_pop), col="grey")
with (ang, text(us_relations ~ church_pop, labels = x, pos = 1, cex=.7) )  
dev.off()

###########################################################################
# Figure 7 (Appendix): Anglican Membership Size/Pop and Relationship w US #
###########################################################################

ang$church_pop_pct <- ang$church_pop_2016 / ang$country_pop_2016
ang <- ang[order(-ang$church_pop_pct),]
x <- c(ang$names[1:3], rep("", 3), ang$names[7:9], rep("", 3) )
y <- c(rep("", 3), ang$names[4], rep("", 8) )
z <- c(rep("", 5), ang$names[6], rep("", 6) )

pdf(file="Fig_7_Ang_Pop_Relationship_Appendix.pdf",width=5,height=4)
plot(ang$church_pop_pct, ang$us_relations, cex.main = .9, ylim = c(.8, 3.3), xlim = c(0, .35), bty="n",
     xlab = "Church population (percentage of national population, 2016)", ylab = "",   yaxt="n",
     main = "African Anglican Church relations with US Church\n(by church population relative to national population)", cex.main=.95
)
axis(2, at = 1:3, labels = c("Relationship\nMaintained", 
                             "Relationship\nStrained",
                             "Relationship\nEnded"), las=2, cex.axis=.7)
with (ang, text(us_relations ~ church_pop_pct, labels = x, pos = 1, cex=.7) )  
with (ang, text(us_relations ~ church_pop_pct, labels = y, pos = 3, cex=.7) )  
with (ang, text(us_relations ~ church_pop_pct, labels = z, pos = 2, cex=.7) )  
dev.off()

##########################################
##########################################
# Figures 3A & 3B from Lutheran Data     #
# Data: Compiled by S. Dreier, 2013-2017 #
##########################################
##########################################

luth_base <- read.csv("luth_data.csv", header=TRUE, stringsAsFactors=F)

#############################################################
# Figure 3(A): Lutheran Membership Size and Numb Companions #
#############################################################

luth <- luth_base[order(-luth_base$numb_elca_companion_2012, -luth_base$church_pop_2012),]
luth$pop <- luth$church_pop_2012/1000000
luth$comp12 <- luth$numb_elca_companion_2012
x <- c(luth$Country[1:7], rep ("", 30))

pdf(file="Fig_3A_Luth_Pop_Compan.pdf",width=5,height=4)
with(luth, plot(jitter(pop), jitter(comp12), cex.main=.9, bty="n",
                xlab="Church population (in millions, 2012)",
                ylab = "Number of ELCA Companion Synods (2012)",
                main = "Number of Companion Synods with\nAfrican Lutheran Churches (by church population)",
                xlim=c(0,5.9),
                ylim=c(-.2,22)
))
abline(lm(luth$comp12 ~ luth$pop), col="grey")
with (luth, text(comp12 ~ pop, labels = x, pos=3, cex=.65) )  
dev.off()

#############################################################
# Figure 3(B): Lutheran Membership Size and US Contribution #
#############################################################

luth <- luth_base[order(-luth_base$church_pop_2012),]
luth$pop <- luth$church_pop_2012/1000000
luth$us_contrib <- luth$elca_companion_usd/1000000
x <- c(luth$Country[1:5], rep("", 32))

pdf(file="Fig_3B_Luth_Pop_USContrib.pdf",width=5,height=4)
par(mai=c(1,1.1,1,.7))
with(luth, plot(jitter(pop), jitter(us_contrib), cex.main=.9, bty="n",
	xlab="Church population (in millions, 2012)",
	ylab = "Investments from ELCA Companion Synods\n(in millions USD, 2008)",
	main = "Companion Synod Contributions\nto African Lutheran Churches (by church population)",
  xlim=c(0,5.9),
  ylim=c(-.3,4)
      ))
abline(lm(luth$us_contrib ~ luth$pop), col="grey")
text(luth$us_contrib ~ luth$pop, labels = x, pos=1, cex=.7) 
dev.off()

##################################################
##################################################
# Figures 4-6 (Appendix): Survey Evidence        #
# Data: Afrobarometer Round 6 (2016)             #
##################################################
##################################################

suppressWarnings(afro6 <- spss.get("merged_r6_data_2016_36countries2.sav", use.value.labels=FALSE))

# Recode for country names
afro6$ctry <- NA
afro6$ctry[afro6$COUNTRY==1] <- "Algeria" #LGBT not asked
afro6$ctry[afro6$COUNTRY==2] <- "Benin" 
afro6$ctry[afro6$COUNTRY==3] <- "Botswana" #Ang, Luth
afro6$ctry[afro6$COUNTRY==4] <- "Burkina Faso"
afro6$ctry[afro6$COUNTRY==5] <- "Burundi" #Ang
afro6$ctry[afro6$COUNTRY==6] <- "Cameroon" #Luth
afro6$ctry[afro6$COUNTRY==7] <- "Cape Verde" 
afro6$ctry[afro6$COUNTRY==8] <- "Cote d'Ivoire"
afro6$ctry[afro6$COUNTRY==9] <- "Egypt" #LGBT not asked 
afro6$ctry[afro6$COUNTRY==10] <- "Gabon" 
afro6$ctry[afro6$COUNTRY==11] <- "Ghana" #Ang, Luth
afro6$ctry[afro6$COUNTRY==12] <- "Guinea" #Ang
afro6$ctry[afro6$COUNTRY==13] <- "Kenya" #Ang, Luth
afro6$ctry[afro6$COUNTRY==14] <- "Lesotho" #Ang
afro6$ctry[afro6$COUNTRY==15] <- "Liberia" #Ang, Luth
afro6$ctry[afro6$COUNTRY==16] <- "Madagascar" #Ang, Luth
afro6$ctry[afro6$COUNTRY==17] <- "Malawi" #Ang, Luth
afro6$ctry[afro6$COUNTRY==18] <- "Mali" 
afro6$ctry[afro6$COUNTRY==19] <- "Mauritius" #Ang
afro6$ctry[afro6$COUNTRY==20] <- "Morocco" 
afro6$ctry[afro6$COUNTRY==21] <- "Mozambique" #Ang, Luth
afro6$ctry[afro6$COUNTRY==22] <- "Namibia" #Ang, Luth
afro6$ctry[afro6$COUNTRY==23] <- "Niger" 
afro6$ctry[afro6$COUNTRY==24] <- "Nigeria" #Ang, Luth
afro6$ctry[afro6$COUNTRY==25] <- "Sao Tome & Principe" 
afro6$ctry[afro6$COUNTRY==26] <- "Senegal"  #Luth
afro6$ctry[afro6$COUNTRY==27] <- "Sierra Leone" #Ang, Luth
afro6$ctry[afro6$COUNTRY==28] <- "South Africa" #Ang, Luth
afro6$ctry[afro6$COUNTRY==29] <- "Sudan" #Ang, LGBT not asked
afro6$ctry[afro6$COUNTRY==30] <- "Swaziland" #Ang
afro6$ctry[afro6$COUNTRY==31] <- "Tanzania" #Ang, Luth
afro6$ctry[afro6$COUNTRY==32] <- "Togo"
afro6$ctry[afro6$COUNTRY==33] <- "Tunisia" 
afro6$ctry[afro6$COUNTRY==34] <- "Uganda" #Ang
afro6$ctry[afro6$COUNTRY==35] <- "Zambia" #Ang, Luth
afro6$ctry[afro6$COUNTRY==36] <- "Zimbabwe" #Ang, Luth

all(table(afro6$ctry) == table(afro6$COUNTRY)) #check recoding accuracy

## Calculate country statistics ##

# LGBT (attitudes)
afro <- subset(afro6, Q89C!=99) #99: not asked in Algeria, Egypt, Sudan
afro$lgbt <- afro$Q89C
afro$lgbt[afro$lgbt==-1 | afro$lgbt==9 | afro$lgbt==99] <- NA

# Religiosity (reported behavior)
afro$religiosity <- afro$Q98B
afro$religiosity[afro$religiosity==-1 | afro$religiosity==9 | afro$religiosity==98] <- NA
afro$religiosity <- afro$religiosity+1
afro$religiosity[afro$religiosity==8] <- 0

# Religious leaders corrupt (attitudes)
afro$rel.corrupt <- afro$Q53I
afro$rel.corrupt[afro$rel.corrupt==-1 | afro$rel.corrupt==9] <- NA

# Trust religious leaders (attitudes)
afro$rel.trust <- afro$Q52L
afro$rel.trust[afro$rel.trust==-1 | afro$rel.trust==9] <- NA

# Contact religious leaders (reported behavior)
afro$rel.contact <- afro$Q24F
afro$rel.contact[afro$rel.contact==-1 | afro$rel.contact==9] <- NA

# Split data by country for country averages
afro_data <- split(afro, afro$ctry)

# Dislike LGBT (percent by country)
dislike_lgbt <- function(dataset){
    round( (table(dataset$lgbt)[1] + table(dataset$lgbt)[2]) / sum(table(dataset$lgbt)) ,3)
    } 

lgbt2 <- lapply(afro_data, dislike_lgbt)
x <- as.numeric(unlist(lgbt2))
y <- as.character(sort(unique(afro$ctry))) # Sort to match lapply order

data <- data.frame(y, x)
data$y <- as.character(data$y)
colnames(data) <- c("Country", "lgbt")

# Religiosity (mean by country)
religiosity_mean <- function(dataset){
  round( mean(na.omit(dataset$religiosity)) ,3)
} 
relig <- lapply(afro_data, religiosity_mean)
relig <- as.numeric(unlist(relig))

# Percent (by country) who say no religious leader is corrupt
corrupt <- function(dataset){
  round( table(dataset$rel.corrupt)[1]  / ( sum(table(dataset$rel.corrupt))) ,3)
} 

rel.corrupt <- lapply(afro_data, corrupt)
rel.corrupt <- as.numeric(unlist(rel.corrupt))

# Percent (by country) who say they trust religious leaders a lot
trust <- function(dataset){
  round( table(dataset$rel.trust)[4]  / ( sum(table(dataset$rel.trust))) ,3)
} 

rel.trust <- lapply(afro_data, trust)
rel.trust <- as.numeric(unlist(rel.trust))

# Percent (by country) who say they consult a  religious leaders often
contact <- function(dataset){
  round( table(dataset$rel.contact)[4]  / ( sum(table(dataset$rel.contact))) ,3)
} 

rel.contact <- lapply(afro_data, contact)
rel.contact <- as.numeric(unlist(rel.contact))

data <- data.frame(data, relig, rel.corrupt, rel.trust, rel.contact) #country means/percentages in one dataset

##################################################
# Figure 4 (Appendix): LGBT Attitudes by country #
##################################################

data0 <- data[rev(order(data$lgbt)),]

pdf(file="Fig_4_LGBT_Ctry_Appendix.pdf",width=5,height=4)
par(mai=c(1.2,1.3,0.7,.8))
plot(data0$lgbt, nrow(data0):1, type="p", xlim =c(.19,1.01), ylim=c(0,33),
     pch=16, yaxt="n", ylab="", bty="n", cex.main=.9, cex.axis=.9,
     main="Public opposition to sexual minorities\n(by country)", xlab="")
axis(2, at=1:33, labels = rev(data0$Country), las=1, cex.axis=.6) 
abline(h=8.5, col="grey")
mtext("% reporting they would dislike having an LGBT neighbor", side=1, line=2.5, cex=.9)
mtext("Source: Afrobarometer Round 6 (2016)", side=1, line=4, cex=.75)
dev.off()

################################################################
# Figure 5 (Appendix): LGBT Attitudes & Religiosity by country #
################################################################

pdf(file="Fig_5_LGBT_Religiosity_Appendix.pdf",width=5,height=4)
plot(data$lgbt, data$relig, type="p", pch=16, bty="n", cex.main=.9, cex=.9, 
     xlab="", ylab="", ylim = c(2.8, 7), yaxt="n",
     main="Public opposition to sexual minorities\nand mean religiosity by country")
axis(2, at=3:7, las=1, cex.axis=.65, 
    labels = c("Once a month",
               "Once a week", 
               "Few times\na week", 
               "Once a day",
               "More than\nonce a day")  )
mtext("% reporting they would dislike having an LGBT neighbor", side=1, line=2.5, cex=.9)
mtext("Frequency of Religious Practice", side=2, line=4.5, cex=.9)
mtext("Source: Afrobarometer Round 6 (2016)", side=1, line=4, cex=.75)
dev.off()

####################################################################
# Figure 6 (Appendix): LGBT Attitudes & Anglican Relationship w US #
####################################################################

ang_merge <- merge(ang_base, data, by.x.y="Country") # Merge Afrobarometer with Anglican data

pdf(file="Fig_6_LGBT_Ang_Relat_Appendix.pdf",width=5,height=4)
plot(ang_merge$lgbt, ang_merge$us_relations, type="p", pch=16, bty="n", cex=.9, cex.main=.9,
     yaxt="n", xlab="", ylab="",
     main="Public opposition to sexual minorities\nand Anglican relationship with US church",
     xlim =c(.29,1.01), ylim=c(.75,3.25)
     )
axis(2, at=1:3, las=1, cex.axis=.65, 
     labels = c("Relationship\nMaintained", 
                "Relationship\nStrained", 
                "Relationship\nEnded")  )
mtext("% reporting they would dislike having an LGBT neighbor", side=1, line=2.5, cex=.9)
mtext("Source: Afrobarometer Round 6 (2016)", side=1, line=4, cex=.75)
dev.off()
